clear
rng default
load IdSetM
load IdSetW
load U_CS
load V_CS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% The core %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ntypes=3;
D_interval=cell(ntypes,ntypes);
for i=1:ntypes-1
    for k=i+1:ntypes
        %Men
        M_ii_ik=unique([IdSetM{i}(:,i+1) IdSetM{i}(:,k+1)], 'rows', 'stable'); 
        M_kk_ki=unique([IdSetM{k}(:,k+1) IdSetM{k}(:,i+1)], 'rows', 'stable');
        [cx1, cx2] = ndgrid(1:size(M_ii_ik,1), 1:size(M_kk_ki,1));
        cx1=cx1(:);
        cx2=cx2(:);
        M_ii_kk_ik_ki=[M_ii_ik(cx1,1) M_kk_ki(cx2,1) M_ii_ik(cx1,2) M_kk_ki(cx2,2)]; 
        
        %Women
        W_ii_ki=unique([IdSetW{i}(:,i+1) IdSetW{i}(:,k+1)], 'rows', 'stable'); 
        W_kk_ik=unique([IdSetW{k}(:,k+1) IdSetW{k}(:,i+1)], 'rows', 'stable');
        [cx1, cx2] = ndgrid(1:size(W_ii_ki,1), 1:size(W_kk_ik,1));
        cx1=cx1(:);
        cx2=cx2(:);
        W_ii_kk_ik_ki=[W_ii_ki(cx1,1) W_kk_ik(cx2,1) W_kk_ik(cx2,2) W_ii_ki(cx1,2)];
        
        %Compose
        %A(h,1)+B(:,1)+A(h,2)+B(:,2)-(A(h,3)+B(:,3))-(A(h,4)+B(:,4))
        BB = W_ii_kk_ik_ki(:,1) + W_ii_kk_ik_ki(:,2) - W_ii_kk_ik_ki(:,3) - W_ii_kk_ik_ki(:,4);
        AA = M_ii_kk_ik_ki(:,1) + M_ii_kk_ik_ki(:,2) - M_ii_kk_ik_ki(:,3) - M_ii_kk_ik_ki(:,4); 
        D_interval{i,k}= [min(BB) max(BB)] + [min(AA) max(AA)];
    end
end

D_CS=zeros(ntypes,ntypes);
for i=1:ntypes-1
    for k=i+1:ntypes
        D_CS(i,k)=U_CS(i,i)+V_CS(i,i)+U_CS(k,k)+V_CS(k,k)-(U_CS(i,k)+V_CS(i,k))-(U_CS(k,i)+V_CS(k,i)); 
    end
end

save('D_interval', 'D_interval')
save('D_CS', 'D_CS')
